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Viable models of modified gravity which satisfy both local as well as cosmological 
tests are investigated. It is demonstrated that some versions of such highly non- 
linear models exhibit multiply de Sitter universe solutions, which often appear in 
pairs, being one of them stable and the other unstable. It is explicitly shown that, 
for some values of the parameters, it is possible to find several de Sitter spaces (as a 
rule, numerically); one of them may serve for the inflationary stage, while the other 
can be used for the description of the dark energy epoch. The numerical evolution of 
the effective equation of state parameter is also presented, showing that these models 
can be considered as natural candidates for the unification of early-time inflation with 
late-time acceleration through dS critical points. Moreover, based on the de Sitter 
solutions, multiply SdS universes are constructed which might also appear at the 
(pre-)inflationary stage. Their thermodynamics are studied and free energies are 
compared. 



PACS numbers: ll.25.-w, 95.36.+X, 98.80.-k 
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I. INTRODUCTION 

Modified gravity (for a review, see e.g. constitutes an interesting dynamical alter- 
native to ACDM cosmology in that it is also able to describe with success the current 
acceleration in the expansion of our Universe, the present dark energy epoch. The specific 



class of modified f{R) gravities (for a review, see e.g. 



2|) has undergone many studies 



which suggest that this family of gravitational alternatives for dark energy is able to pass 



the constringent solar system tests. The investigation of cosmic acceleration as we^ 



1 as the 



study of the cosmological properties of f{R) models has been done in Refs. [ll,l2|,|3|,y,l5|, l6|. 
The possibility of a natural unification of early-time inflation with late-time acceleration 
becomes a realistic and quite natural possibility in such models, as is demonstrated e.g. 
in Ref. fl. 

Recently, the importance of modified gravity models of this kind has been reassessed 
with the appearance of the so-called 'viable' f{R) models 3, Isl, Is, 10, 11|. Those are 
theories which satisfy both the cosmological as well as the local gravity constraints, which 
had caused in the past a number of serious problems to some of the first-generation 
theories, that had to be considered now as only approximate descriptions from more 
realistic theories. The final aim of all those phenomenological models is to describe a 
segment as large as possible of the entire history of our universe, as well as to recover all 
local predictions of Einstein's gravity that have been already verified experimentally, with 
very good accuracy, at the solar system scale. It is remarkable that, as was demonstrated 
in Refs. , some of these realistic models lead to a natural unification of the early-time 

inflation epoch with the late-time acceleration stage. 

Let us recall that, in general (see e.g. , for a review), the total action for the 

modified /(-R) gravitational models can be written as 



S = ^j d'x^ [R - F{R)] + S^m) 



Here F{R) is a suitable function of the scalar curvature R, which defines the modified 
gravitational part of the model. The general equation of motion in f{R) = R — F{R) 
gravity with matter is given by 



lg,uf{R) - R,uf\R) - g,uUf\R) + V,VJ\R) = -yT(„,)^. 



(2) 
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where T{rn)^v is the matter energy- momentum tensor and f'{R) is the derivative of f{R) 
with respect to its argument R. For a generic f{R) model is not easy to find exact static 
solutions. However, if one impose some restrictions, one can proceed along the following 
lines. First, we may require the existence of solutions with constant scalar curvature 
R = Rq, and we arrive at 

f'{Ro)Ri_iu = — Y~giJLv (3) 

Taking the trace, we have the condition 

2f{Ro) = RonRo) (4) 

and this means that the solutions are Einstein's spaces, namely they have to satisfy the 
equation 

-'^f^!^ 2/'(i?o) A '' 

Rq being a solution of Eq. (jll). This gives rise to an effective cosmological constant, 
namely 



^eff = 7777777T = "T • (6) 



f{Ro) _ Rq 
2f'{Ro) ~ 4 

The purpose of our work will be to study the appearance of multiply de Sitter space 
solutions in several realistic models of modified gravity. The occurrence of multiply de 
Sitter solutions plays a fundamental role in modified gravity because it permits to describe 
the infiation stage as well as current ACDM cosmology in terms of this theory alone 
without any need for either fine-tuning of a cosmological constant nor of introducing 
extra scalar fields. In summary, this is a minimal and at the same time very powerful 
approach which circumvents some of the hardest problems of present day physics. One 
may argue that these theories are equivalent to introducing extra scalar fields, but this 
equivalence has been proven to hold at the classical level only, not at the quantum one. 
In addition, the cosmological interpretation of modified gravity solutions is different from 
that of scalar field cosmology. 

The paper is organized as follows. In the next section we discuss a viable modified 
gravity model [3] in an attempt to study de Sitter solutions there. It is shown that, 
for some values of the parameters, it is possible to find several de Sitter spaces (as a 
rule, numerically, for small values of the curvature exponent); one of them may serve for 
the infiationary stage, while the other one can be used for the description of the dark 
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energy stage. The evolution of the effective equation of state parameter is investigated 
numerically. Sect. 3 is devoted to the analysis of the same question in a slightly generalized 
model which is known to describe the unification of the early-time inflation with the 
late-time acceleration epochs j^. The same numerical investigation is carried out, and 
a multiply de Sitter universe solution is constructed. In section 4, the corresponding 
em is investigated for a viable model of tangential modified gravity, proposed in 



prob 

Ref . [Ul] . Having in mind the possibility to use a SdS universe for the description of cosmic 
acceleration, those solutions for modified gravity are investigated in Sect. 5. Comparison 
of the free energies of the de Sitter and the SdS solutions, for the viable model of the 
second section, is then made. Finally, some outlook is given, together with the conclusions, 
at the end of the paper. 



II. DYNAMICAL SYSTEM APPROACH AND DE SITTER SOLUTIONS IN 

REALISTIC MODIFIED GRAVITY 



In this section we will study de Sitter solutions in realistic modified gravity models 
using a dynamical system approach. We shall start from the following form for the initial 
action 



S 



1 



,R-F{A) 



(7) 



where A is some function of the geometrical invariants. In order to make use of the 
method of dynamical systems, we work in the metric corresponding to a spatially-fiat 
FRW universe, namely 



diag{-n{tf, a{tf, a{tf, a{tf). 



(8) 



From here, the FRW equation can be written in the following way (for more details see 



Ref. 



12|) 
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Z^H"^ = Pf, 



where 



PF = F + F'An - SHF'An - F 



,dAn 
'df 



F"A.A. 



(9) 



(10) 



The "dot" over the symbol means derivative with respect to time t, while A^ and An 
represent the derivative of A with respect to n and h respectively. In what follows we will 
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mainly concentrate on the simplest case A = R. 

Q 

Let us consider the following choice for the function F [7|], which represents a very 
interesting subclass of viable modified gravitational models 



F{R) 



2«:2 C2{4r + 1 



where the constant fi has dimension of mass, while ci, C2, C3 are some positive dimen- 



sionless constants. Note that here R = \ 



a 9:lh -i- _ 

a a n 



Rn = — 2i? 



-12{H + 2H^ 



Rn = —6H. After making variation, it is chosen n = 1. Introducing a nonzero constant 
C3 is here equivalent to introducing a shift in the effective cosmological constant, and for 
this reason we will assume C3 = in our further computations. Finally, we can write 
Eq. ([9]) for the function in (|TT|) as 



2\kci ( R 



+ 



(12) 



Z&jHii+AH'^H) kci 



{k 



, k-2 / \ 2fc-2 

f -C2(k + 1) 3 



(13) 



where R = 6{H + 2H^). This equation can be rewritten as a dynamical system, namely 

H = C, 
C = F,{H,C). 

It is easy to see that the critical points of this system are the de Sitter points {H = 0, H = 
0). To investigate the nature of these points we need to determine them explicitly, a non- 
trivial problem in the general case. It is however easy to obtain from (1121) the equation 
satisfied by the critical points Hq: 



7' ' 



Cl 



6K 



2 fcci / 

^2 I ,^2 



k-1 



C2 



12^2 ^ 



(l2Hl 



2 • 



The same result can be obtained directly starting from equation 
Section. 



(14) 

of the previous 



^2 

It is convenient to introduce the notation xq = 12^, for further investigation of this 



equation. So, finally, we have 



^2 2fc+l 

C2X0 



- 2ciC2x'o' + 2c2xl+' + 2ci(-/e - 1)^^ + = 0. 



(15) 
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First of all we find one (trivial) root of this equation xq = (it corresponds to Ho = 0), 
what allows us to reduce the order of the equation. (Note also that, if C3 7^ 0, Eq. ( |T5l) 
takes a more complicate form and then it does not have the trivial root xq = 0.) But 
nevertheless the equation still is of 2A;-order and in the interesting case [k > 2); this is 
too high and the roots cannot be found algebraically. We need a specific discussion of 
this problem. It is clear that if there are too many (say 10 or 20) de Sitter points in 
the theory, it looks like a classical analogue (at a reduced scale, of course) ISj of the 
string landscape vacuum structure in which case it will be far from trivial to obtain the 
standard cosmology. Nevertheless, something can be done even in this case, by comparing 
the energies of the corresponding de Sitter solutions, which should in fact differ, as is 
discussed in Sect. V. Using Descartes rule of signs, we find that Eq. ffT^ can have 2 or 
positive roots. More detailed information can be obtained by using Sturm's theorem. 
Unfortunately it is not possible to apply it in the general case (for arbitrary k), so we 
need to investigate separately each of the different cases, 3 ^ A; ^ 10, and see which one 
is more interesting from a cosmological viewpoint. We have therefore undertaken here a 
systematic analysis of all possible cases in this region of values of k. 

Computations are rather involved (specially for larger k), but the final results are not 
so difficult to describe analytically. We have found, for any of the values of k considered, 
in the range above, that the number of roots depends on the parameter = -khj only, 
and that there is an al such that, for < ak < al, Eq. ( |T5i) has no positive roots, 
while for ak > al, Eq. (fT5l) has two positive roots. We here now enumerate our results 



systematically: a*^ = + ^ 2.43; al = ^ - ^33 ^ 2.01; a*^ = + 

~ 1 4Q- r^,* — 3^17-6977 _ 2-8329-32183 . ^ i no. „* _ 7^17681 , 7^107 . /77T ~ 

352 V J ~ l.'iJ, Wg — 2245 22'l53 V ~ l.UO, LKy — 2-33513 2-32513 V^J- ~ 

0.68; al = ^ + ^ 0.43; a; = ^^^^09 + aiipss^ ^ o.26; a*,, = 

5i°4oi843202307 ^ 5^^17^^9064429^ ^ q ^ote also that, near the critical point a*,, 
both solution are very close to each other when they exist (a real value of ak > al), and 
they are complex conjugate to each other for < afc < al- This means that, in a situation 
where these two roots differ substantially, it must necessarily be ak ^ al- 

The dS-points described above can be used for the construction of inflationary or late- 
time acceleration behavior (depending on the value of the scale factor /i^). To this aim, 
the corresponding dS point must be unstable in the inflationary stage but can be either 
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unstable or stable for the late-time acceleration one. However, this is in fact not a strictly 
necessary condition, since even for stable dS inflation, the exit from it can be achieved 
by a coupling with matter, through the effect of small non-local term or by some other 
mechanism. Unfortunately, an all-round investigation of stability of the dS point turns 
out to be very difficult in the general case. For the recent analysis of critical points in more 
general modified gravity theories depending on all geometrical invariant see Ref. [l^ . 

Here we can carry out the analysis on the stability of de Sitter points for the model we 
are dealing with in this section, only for some chosen numerical set of parameters. Such 
computations show that one of the dS points is very likely to be always unstable (with a 
smaller value of H) and that another one is very likely to be always stable (with larger 
value of H). This means that there is actually an easy way for the inflationary stage 
construction (if the initial conditions lie sufficiently close to the unstable dS point, see 
Fig. [T]). However, only within the model under consideration it might be a problem with 
late-time acceleration, originated from the stable point. The only possibility for late time 
acceleration is to choose very special parameters, in which the stable point is situated very 
far away from the unstable one. Of course, the problem disappears when one describes 
only late-time acceleration within such a model, or when one takes into account other 
terms like the local ones. Note however that specific values of the parameters in which 
both points are unstable can be chosen too, which give rise to a sound theory which 
provides a unified description of the inflationary and late-time acceleration epochs. 

Let us now consider specific numerical results. First, note that a full numerical inves- 
tigation starting from the inflation (or from our dS point) and extending to the late-time 
acceleration epoch and going through a FRW-like stage is not feasible because of the 
presence of numerical instabilities that prevent such possibility. For this reason, we will 
investigate the regions near the equilibrium points only, and especially near the unstable 
one. In Fig. [1] the trajectories in the phase plane {H — H) are presented for the set of 
parameters: fi = l,k = A, Ci = 1.12 ■ 10"'^, C2 = 6.21 ■ 10"^. Such small values of Ci and C2 
are motivated in 7|, but the picture is typical in any case. The position of the dS point 
is specified. We may see one stable point (with bigger value of H) and another, unstable 
one. The stars in the graph correspond to the beginning of the trajectories. We see clearly 
that trajectories can pass sufficiently close to the dS point, and that they escape from 
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FIG. 1: Phase plane (H — C = H) near the equihbrium points. Stars denote the beginning of 
the trajectories. Dashed and dot dashed hnes lead to bottom. 




t 



FIG. 2: Evolution of Weff near an unstable point for the dashed and dot-dashed lines of Fig. [H 
The time origin corresponds to the stars from the previous graph. 

it, lead by the two attractors. One of them (the upper one) leads to a singularity, while 
another one (heading towards the bottom), leads to a solution going somewhere in the 
FRW-like region (this is difficult to analyze exactly, due to numerical instabilities). In 
Fig. [2] the evolution of Wejf = — 1 — is depicted for two trajectories: the upper one 
corresponds to the dashed line of Fig. [1], the bottom one to the dot-dashed line of Fig. [T] 
(both of them lead to the bottom in Fig. [1]). We see that there actually are trajectories 
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which start from phantom-hke solutions and tend to a normal evolution (when this is 
possible in the model). Of course, these are qualitative considerations only, owing to the 
fact that no matter is taken into account. Most probably, inclusion of matter will change 
the situation. However, it is just remarkable that both the early-time inflation epoch as 
well as the late-time accelerating one could be obtained in a unified and natural way in 
such a model, owing to the presence of several dS points. 



III. A DIFFERENT EXAMPLE OF f{R) VIABLE MODEL 

A slightly different choice for the function F{R) is motivated by the realistic and 
viable model which was proposed for the unification of early-time inflation and late-time 
acceleration in Ref. j^, namely 

Here we imply that /o, /i and Rq is positive, what seems reasonable. Now, our FRW-like 
equation Qi reads 



SeiHH+,H^H)foMR-R..r~^ _ ^^^^ _ ^^^^^ _ ^ y 

{fo+fi{(R~Ro)"'+R^}) J'- 0/ \ jJ>-\ 



(17) 



where R = 6{H + 2H'^). Rewriting this equation under the form of a dynamical system 
( IT3l) we find, by introducing the new variable x = I^Hq — Rq, that the equilibrium points 
of this system (which are dS points indeed) are obtained from the following equation 

[2i?o/i(/o + fiK) - 8^^f\R^ + 2t,\m - 2)/o] + 2t^R^Umx^-^+ (18) 

(/o + hKfx + i?o(/o + hK){k + fiK - 4^2/?-"!) = 0. 

It is clear that, in the general case (even for fixed m > 2), solving this equation will not 
be an easy thing. We can make use of the following hint. Let us require that x = be a 
solution of Eq. ( ITSl) . This is possible only if there is some relation among the constants 
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of our theory. As we can easily see from (fTSj) . it must be that /o = -R™~^(4fi:^ — /i-Ro)- 
Moreover, since /o > 0, we need that /i-Rq < 4k^. Substituting this condition into ffTSl) 
we find an essential simplification of the former equation, namely 

(19) 

where the root x = is already excluded. We now discuss the physical meaning of finding 
a solution x = 0. As a rule this will imply that Rq is the value of the scalar curvature at 
present time, but in principle -Rq is just a parameter of the theory, which could be given 
any likely value. We can, for instance, fix Rq to have the dS point value which is situated 
exactly at present time, to, or at to + 100 years. And there is also the possibility to explain 
late time acceleration if this point is stable. Now let us consider other possible roots of 
Eq. ffTUl) . Using Sturm's theorem as in the previous section, we find several positive roots 
(which correspond to dS points in the past) and a number of negative roots (needless to 
say, only values which are x > —Rq have a physical meaning), corresponding to dS points 
in the future. Unfortunately, computations are harder than in the previous case and our 
calculations have been performed for three specific values of the parameter m: 3, 5 and 
7, only. In any case, the results are very similar for the three situations, what gives us 
hope that for larger values of m result will be also similar. We may predict the upper 
limit of the dS points in the future only, which is 2, but in any case this point is not much 
interesting. 

Concerning the dS points in the past, the picture turns out to be much like the one in 
the previous case: the number of dS points depends only on a dimensionless parameter 
for all investigated values of m. That is, P = And there are two values, and 

P^^, which are different for each m, so for < /3 < and /5 > there are two dS 
points and for [3^, < [3 < there is no dS point. We have not been able to obtain exact 
analytical expressions for and and thus give here numerical results only. For m = 3 

^ 0.153, ^ 3.12 (as positive roots of A[3^ + 336/5=^ - 843/32 - 815/3 + 144 = 0). For 
m = 5/3* ^ 0.163, /3,* ^ 2.01 (as positive roots of 1318032/3^40509072/3^-208593144/?^+ 
472402800/3^-637548615/3^+422702939/3=^-90 8 9 631/32 - 69 1 90 9 83/3+1230 1 875 = 0). For 
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m = 7 ^ 0.215, ^ 1.68 (as positive roots of 2^5^24337/3^^ + 2^5 ■ 181772596417/3^° - 
2^7 ■ 17 ■ 29 ■ 641 ■ 15074567/3^ + 2^19 ■ 29 ■ 173 • 6421 ■ 100207/3^ - 2^7 ■ 153140941633579/3^ + 
2^7 ■ 227057313819467/3*^ - 5 ■ 7 ■ 195531002131489/3^ + 2^47 ■ 1331611 ■ 21800759/3^ - 2 ■ 5 ■ 
53 ■ 2797 ■ 2102824303/3^ + 2^7 ■ 29 ■ 728096696819/3^ - 5-^23 ■ 18382295597/3 + 2^^5^7^ = 0). 

As we can see, the interval where there are no dS points shrinks when the value of the 
parameter m is increased. If we assume that for larger m the number of dS points depends 
on the parameter (3 only — as clearly happens in the investigated cases — we do find /3* and 
for any value of m. Such investigations show that /3* and /3** slowly change when m 
increases and that the interval without dS points is present even for very big value of the 
parameter m. For example, for m = 101 we have /3* ~ 0.480 and /3** ~ 1.06. A simple 
numerical investigation of stability for the existing dS points shows that one of them (the 
one with the smaller value of H) is always unstable, but the other one (with the larger 
value of H) can either be stable or unstable. From our numerical results, the appearance 
of unstable dS points pairs relative to the interval < /3 < /3* seems most probable. Note 
also that the condition /o > means that it actually must be /3 > 0.5 and that in this 
case it is most likely the one stable and one unstable point show up. 

Now let us consider the evolution equation 0171) from a different point of view. As we 
already noted, Rq is a parameter of the theory which corresponds to the value of the scalar 
curvature at some epoch. This means that during the normal evolution of our universe, 
from large to little (or zero) curvature, it becomes R = Rq aX some moment. But, as 
we can see from Eq. (ITTl) . this means that that at this point the coefficient of the higher 
derivative term {H) is equal to zero. This is a well-known mathematical problem, which 
needs special investigation. As we know from mathematics there are two possibilities: 
the solution of the perturbed equation (with a higher derivative term) may tend to the 
solution of the degenerate equation (without higher derivatives), and then the coefficient 
of the higher derivative term may tend to zero or not tend to the solution of the degenerate 
equation. A special investigation of this problem shows that solutions of Eq. (fT7|) tend to 
solutions of the degenerate equation, which is 

= T^^^. (20) 



K 



ifo + fiKY 



when R = Rq + O but not so when R = Rq — O. This means that the point Rq is reachable 
during the evolution from large R to zero, but it is not the final point of the evolution. 
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because there is an instability in the future. On the other hand, we have R = Rq ai this 
point, and therefore 3if| = i?o/4. Substituting this into ( !20l) we find a relation among 
the parameters of the theory which were introduced before, /o = -R™~^(4/t^ — fiRo). This 
means that, generically, if one wants the point Rq to be reachable during evolution, a 
relation of this sort among the parameters of the theory must be fulfilled. That is, we 
have two independent parameters only (or even just one if we consider that Rq is strictly 
related to our epoch). Note also that it is impossible to use standard numerical methods 
near this point for solving Eq. (fT7|) because all available methods need solving the equation 
with respect to the highest derivative, and it turns out that, near this point, the numerical 
solution is unstable. 



0.02 
0.00 

c 

-0.02 
-0.04 

0.43 0.48 0.53 0.58 

H 

FIG. 3: Phase plane {H—C = H) near equilibrium points. Stars denote beginning of trajectories. 
Dashed and dot dashed lines lead to bottom. 

Now let us consider numerical results near the stable points. Since the value of Rq is 
very small in our epoch (~ 10"^^ cm~^), it is impossible to use real values in computations. 
Thus, we have produced a qualitative analysis using the following set of parameters: 
Rq = 1, /i = 1, = 1, m = 7. As we can see, (3 = 2, and there must be two dS 
points corresponding to positive x. The evolution lines near these points are represented 
in Fig. [3l where we can trace one stable and one unstable point, stars denoting the 
beginning of trajectories, as in the previous case. Actually, there are three dS points 
for these chosen parameters. One of them, H = 0.288, which corresponds to x = 0, is 
degenerate and situated out of the graph. The time evolution of We/f is represented in 
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-0.95 



-1.10 




FIG. 4: Evolution of w^jf near unstable point for dashed and dot-dashed lines of Fig(TJ Zero in 
time corresponds to the stars from the previous graph. 

Fig. m by the dashed and dot-dashed Unes. The picture is very similar to the one in the 
previous section, and the same comments also apply here. 

Thus, we have here shown that the unification of early-time inflation with late-time 
acceleration is in principle possible, and even very likely, due to the appearance of several 
dS points in the evolution of the universe. For a more realistic study, the presence of 
(minimally and non-minimally coupled) matter should be taken into account with care. 



IV. DE SITTER UNIVERSE FROM TANGENTIAL MODIFIED GRAVITY 



As a third example, now we discuss a model proposed in ll|]. Such model is deflned 
by means of the function 



f{R) =R- F{R) = R-a 



tanh ( — — z ^ I + tanh ° 



f211 



2 / V 2 

where a, b and Rq are arbitrary parameters. One immediately sees that -F(O) = 0, as 
required and, moreover, that 

'bRo' 



lim F{R) = 2Aeff = a 

R— »oo 



1 + tanh 



(22) 
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If i? ^ -Ro in the present universe then Acs plays the role of the effective cosmological 
constant. We also observe that the derivative 

f'{R) = 1 f ^ (23) 

2cosh^(^(^' 

has a minimum when R = Rr,, which reads 



f'iRo) = 1 - y • (24) 

In order to avoid antigravity one needs to require 

< f'iR) < f'iRo) = 1 - y • (25) 

The model given by Eq. (1211) is able to describe late acceleration, since, in general, de 
Sitter critical points exist. They are the solutions of (jlj). Then, we set 

KiR) = 2f{R) - Rf\R) , (26) 

and study (numerically) the zeros of this transcendental function (see Fig. [5]). We see 
that for suitable choices of the parameters there are one or two de Sitter critical points. 




FIG. 5: K{R) for Rq = l,a = 1.64 and different choices of the parameter b: b = 1, b = a (picture 
on the left); b = —a, b = —10a (picture on the right). 

Now, we impose Rq to be a de Sitter critical point for the model in Eq. fl?T|) . This 
means that (j26l) has to be satisfied for R = Rq and such a condition fixes one parameter. 
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say a. Then, we get 



2Rn 



bRo-4:tanh{bRo/2) ' 
and introducing, for convenience, the dimensionless variables 

R 



X 



R(] 



bo = Rob , 



we obtain 



/'(i?o) = 1 



60 - 4tanh(6o/2) 



A 



eff 



1 + tanh(6o/2) 
60 -4tanh(6o/2) ' 



(27) 



(2^ 



(29) 



^ ^ (1 - x) tanh(&o/2) + 4tanh(6o(3: - l)/2) + 6ox tanh'(6o(x - l)/2) 

^ ' 60 -4tanh(6oa;/2) ' ^ ' 

In order to have Acs > and, at the same time, to avoid antigravity, the parameter bo 
has to be negative and in the range bo < 4tanh(6o/2) ~ —3.83. By varying bo in that 
range, the value of A^s can acquire any desired value. Moreover the model has always two 
de Sitter critical points with constant curvatures Ri and R2. One can see that R2 = Ro, 
while Ro/2 < Ri < Rq. In Fig. 1^ we have plotted K{R) for some values of bo. Thus, our 





0.2 0.4 0.6 0.8 1 



FIG. 6: K{R), with Ro as a fixed de Sitter solution, for diff'erent choices of the parameter 
bo = —3.84, bo = —3.85 (picture on the left); 60 = —10, bo = —20 (picture on the right). 



study shows the existence of multiply-de Sitter universes in more complicated, tangential 
models. One can prove that other viable F(R) gravities introduced in ref. [ll| naturally 



lead to de Sitter solutions too, as demonstrated in Ref. 15 1 
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V. BLACK HOLE SOLUTIONS AND RELATED THERMODYNAMICAL 

QUANTITIES 

In this section we will discuss spherically symmetric exact solutions for the modified 
gravities above. We mainly concentrate on multiply Black Hole solutions and on their 
related entropies and free energies. 

We have seen in the introduction that, requiring constant curvature solutions, one may 
consider the simplified form of the equations of motion given in ([5]). As a result, f{R) 
modified gravity models admit the following class of general static neutral black hole 
solutions — in four dimensions and with a non vanishing cosmological constant. In order 
to describe them, recall the metric 

ds^ = -V{r)de +y^ + r^dEl , V{r) > , (31) 

where k = 0, ±1 and the possible horizon manifolds are Si = S"^, the two dimensional 
sphere, So = T^, the two dimensional torus, and S_i = H^/T, the two dimensional 
compact Riemann surface. The scalar curvature for the ansatz (l3Tl) reads 

R = -l.[r'^V" + ArV' + 2V -2k] (32) 

and this means that every constant curvature solution with R = Rq has to satisfy the 
equation 

r V" + ArV + 2V-2k = -r^Ro . (33) 

The general solution of this differential equation depends on two integration constants, b 
and c, and reads 

V{r) = \ + k-'--^r\ (34) 
The Ricci tensor and the scalar curvature have to satisfy the equations 



Rq 

Rfjiv = Qiiv , R = Ro, (35) 



where Rq is implicitly given by (jl]), that is 



2/(i?o) 
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Now, it is easy to verify that Eqs. ( l35l) are satisfied only if 6 = 0, while c is an arbitrary 
parameter, which is usually assumed to be non-negative and related to the mass M of the 
black hole by c = 2MG. The special case c = is also admissible. Then, we have 

V(r) = k---—r\ c>0, V(r)>0. (37) 
r 12 

In principle there are physical solutions of the latter equation for k = and k = —1, 
which give rise, respectively, to a torus topology and a hyperbolic topology for the horizon 
manifold (the so called topological black holes 16| ) . Here we are mainly interested in the 
usual spherical symmetric horizons and so we only consider in detail the case k = 1. 

As it is well known, in the special case k = l,c = 0, V{r) in ( |371) is always positive 
when Rq < 0, and this corresponds to the Anti de Sitter (AdS) solution. On the contrary, 
when Rq > 0, V{r) is positive, for r < 2^/3jRo, and this corresponds to the de Sitter 
solution. If k = 1 and c = 2MG > 0, one has black hole solutions but only if 

c^Ro - 7r<0 =^ a = - A/G < 1 • (38) 
9 2 

We see that (l38l) is always satisfied if Rq < 0. This corresponds to the Schwarzschild- 
Anti-de Sitter (SAdS) black hole. In this case r > th, th being the positive root of 
V{r) = (horizon radius). If -Rq > 0, there are solutions only if a < 1. In this case 

^ ^c, fH and rc being the positive roots of V{r) = (horizon and cosmological 
radius respectively). The extremal case a = 1 is also admissible (Nariai solution), but 
the thermodynamics of such a black hole have to be discussed separately 

At this point, we provide a brief discussion regarding the thermo dynamical properties 
of the above black hole solutions. If one make use of the Noether charge method for 
evaluating the entropy associated with the black hole solutions with constant curvature 
jisj l in modified f{R) gravity models, one has [3] 

S = ^fiRn) . (39) 

where the factor f{R) has to be evaluated on the horizon, with area Ah- As a consequence, 
one obtains a modification of the "Area Law". Several examples have been discussed in 
jigl ]. In the case of a constant curvature solution one has simply Rh = Ro, Ro being 
the solution of Eq. (jlj), which has been investigated in previous Sections. With regard to 
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this, we have found sufficient conditions to have two de Sitter solutions. Thus, we may 
investigate their thermodynamical behavior evaluating the associated free energy. 
The free energy T is related to the canonical partition function Z by 

^=-!^. (40) 

On the other hand, a semiclassical approximation gives 

Z ~ e~^^ , (41) 

where Ie is the Euclidean classical action associated with the de Sitter solution. A direct 
calculation leads to 



2471 f{R 







and from (jl]) and (IH^ it directly follows that 



Ie = -Sh^J'=-^, (43) 

Ph 

Sh and Ph being, respectively, the entropy and the inverse temperature of the black hole. 
For a generic Schwarzschild-de Sitter (SdS) solution, Eq. fHHj) reads 

Grio 

The temperature is related to the horizon radius r^^ by 

2MG ruRo 



1 _ \V'{rH)\ 1 



Ph 47r 4:71 
th being a positive solution of the algebraic equation 



r]j 6 



(45) 



As a result, we finally have 



. 12 24MG , , 

4-— rH + ^— = 0. (46) 



T = -\th- 3MG| 4^ • (47) 
Gilo 



In the pure de Sitter case M = 0, rj^ = 2^3/i?o, and so one has 



^^s = -2V3^,. {4i 



GRq 
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For SdS one needs to consider separately the two admissible cases a = 1 and a < 1. 
Here we discuss the second one only, that is, the proper SdS solution where a = 
(3/2)A'/G-\/Ro < 1. In this case, Eq. ( l46l) has one negative root and two distinct positive 
roots, th, the event horizon and rc, the cosmological horizon, with th < rc- The positive 
roots can be written in the form 



rc 



where 



and 



7c= 2 



47c 
47h 

1 



1 



< 7c < 1 



0<lH < 



2 ' 



a 



,1/3 



a 



1/3 



7h 



7c -7c 



In this way, we get 



c 



H 



|27c - a\ 
|27g - «| 



dS ■ 



dS ■ 



(49) 
(50) 

(51) 

(52) 

(53) 
(54) 



Tds being the expression 0481) . It has to be noted that, in the latter equations, the 
factors in front of Tds are always smaller than 1 and thus the energies Tc and Th are 
always smaller than Tds, independently of the value of the mass. Moreover, the factor 
1 27(7 — Oi\l\f?) is a monotone decreasing function of a, which is equal to 1 for a = and 
to for a = 1, while the factor |27// — a\l \f?) is exactly equal to for a = 0, 1 and quite 
small otherwise. Its maximum value, a little bit more than 1/10, is reached for a ~ 0.8. 
This means that the minimum value for the free energy Tu of a SdS black hole is obtained 
for M ~ 0.53/(G'/R;;). 

Now, we shall study some models and explicitly compute the corresponding free energy. 
First of all, as a trivial example, we consider the h^CDM model described by f{R) = 
R — 2A. This has a stable de Sitter solution, with Rq = 4A. Then, we immediately have 

%/3 



dS 



2GVA 



(55) 
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FIG. 7: Tc and Tu SdS free energies for the ACDM model. 

The corresponding SdS cosmological and black-hole free energies are plotted in Fig. [71 as 
a function of the mass M (in energy units G-\/A). 

As a second example we consider the viable model in fill I) , with = 4, ci = 1.12- 10"^, 
C2 = 6.21 ■ 10"^, C3 = 0, as discussed in Sect. II. With this choice of parameters, Eq. (j4j) 
has three real solutions for i?o = (Minkowski), R\ ~ 15.33/i^ and R2 — 34.92/i^. The 
last one (i?2) is stable, while the other two are unstable. For the free energies of the two 
de Sitter solutions we get 

_ _^ /™ ^ _ 0^ ^(2) - _2v/3 ^ ~ (56) 



We see that 

^(2) 

T^l < :F^^^ , -f ~ 3.636 . (57) 

independently of the parameter // and so the stable de Sitter solution with curvature R2 is 



favorite, from the energetic point of view, with respect to the one with smaller curvature 

The SdS free energies are plotted in Fig. [HI as functions of the mass M (in energy 
units 1/fiG), for < M < 2/ {3Gy/Ri). However, it has to be stressed that, for the SdS 
corresponding to the de Sitter solution with curvature R2, the mass needs to be restricted 
to < M < 2/ (3Gv^^) < 2/(3G'v^^) and thus, for a fixed mass M, we can compare the 
corresponding free energies in this range only. For both the cosmological and black-hole 
horizons there are critical values of the black hole mass M, say Mc < Mh, for which 
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J^Q^'fj = ^c^'h- Then, we can distinguish three different regions: 

(a) M < Mc'- in this first case J-'q^h < ^'c'h'i 

(b) Mc < M < Mr: in this second case J^^^^ < J^^^\ but J^^^^ > J^^j}^; 

(c) finally, for 2/{3G^/R^) > M > MH,we have J^^Jii > ^c^h- 

As we have seen from general considerations, the free energies corresponding to 
the cosmological horizons are monotonous functions of the mass, while the ones corre- 
sponding to the black-hole horizons are convex functions, which reach the minimum at 
M ~ 0.53/ {Gy/Ro). For our two cases, Ri and R2, we get 

0.53 0.135 0.008 

0.53 0.089 <2) 0.032 , , 

' G^fR2 ^G ' /iG • ^ ^ 

Then we see that, for the model in ( ITTi) we are dealing with, the configuration with the 
minimum free energy corresponds to SdS with mass M = M2 ~ 0.089 /{fiG) and curvature 
i? = i?2 ~ /34.92/i^. The associated free energy is given in fl59p . We also observe that 
this is smaller than the free energy of the pure de Sitter configuration as it is clear from 



Thus, it is demonstrated here that multiple de Sitter solutions can appear also under 
the form of SdS solutions. In other words, the number of multiply solutions becomes 
significantly bigger. With the appearance of both dS and SdS universe solutions one can 
suggest various scenarios for the universe evolution. For instance, we can conjecture that 
a (pre-) inflationary universe is described by some SdS spacetime. As time proceeds, this 
universe decays and enters into the well known radiation/matter dominance phase. In 
its further evolution, the universe transits to the dS (or almost dS) era, by stability and 
least-energy principle considerations. The future universe may again appear as an SdS 
spacetime. 



VI. DISCUSSION 



In summary, we have investigated in this paper several viable models of modified 
gravity which satisfy both the constraints of local as well as cosmological tests. By means 
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FIG. 8: SdS free energies for the two de Sitter critical points of the model (jlip : cosmological 
free energies on the left and black- hole free energies on the right. 

of a numerical study, it is demonstrated that some versions of highly non-linear models 
exhibit multiply de Sitter universe solutions which often appear in pairs, being one of 
them stable and the other unstable. The numerical evolution of the effective equation of 
state parameter is presented too. As a result, these models can be considered as natural 
candidates for the unification of early-time inflation with late-time acceleration through 
dS critical points. Moreover, based on the de Sitter solutions, multiply SdS solutions can 
also be constructed. Further, we have investigated the thermodynamic properties of these 
SdS universes: their corresponding entropies and free energies have been calculated and 
compared. SdS universe might also appear at the (pre-)inflationary stage. 

Owing to the highly non-linear structure of the theories under discussion, the dS uni- 
verses had to be constructed with numerical tools mainly. Moreover, in order to simplify 
the problem, at this first stage of the investigation we did not to consider matter contribu- 
tions. It is clear that, in the next step, we must necessarily include matter and reconsider 
the problem in its presence. This has the potential to lead to a sufficiently realistic quan- 
titative description of the universe expansion history, in which modified gravity would 
be responsible for both acceleration stages: the infiation epoch and the dark energy one. 
This quantitative analysis will be presented elsewhere. 
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